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Abstract We study DNA self-assembly and DNA computation using 
a coarse-grained DNA model within the directional dynamic bonding 
framework [C. Svaneborg, Comp. Phys. Comm. 183, 1793 (2012)]. In our 
model, a single nucleotide or domain is represented by a single interac- 
tion site. Complementary sites can reversibly hybridize and dehybridize 
during a simulation. This bond dynamics induces a dynamics of the 
angular and dihedral bonds, that model the collective effects of chem- 
ical structure on the hybridization dynamics. We use the DNA model 
to perform simulations of the self-assembly kinetics of DNA tetrahedra, 
an icosahedron, as well as strand displacement operations used in DNA 
computation. 

1 Introduction 

Sequence specific hybridization of DNA single strands makes DNA molecules a 
flexible programmable building block. By choosing the right sequences, DNA self- 
assembly behavior can be programmed to produce well defined nano-structures. 
In the pioneering work of Seeman et al., branched DNA constructs have been uti- 
lized to self- assemble into a variety of structures [32 8 36 35J. With DNA origamis 
Rothemund invented a way to fold long DNA single strands into well defined pla- 
nar structures by adding a large number of short stabilizing oligomer strands [26J. 
Later it was demonstrated how to let the planar origamis self-assemble into 3D 
nano-structures such as a box [2]. Ever since the pioneering work of Adlemann 
in 1994 [lj, DNA has also been recognized as a massively parallel, versatile, and 
inexpensive computing substrate. In order for such substrate to be of practical 
interest, however, it is desirable that the computational framework is scalable 
and that individual computational elements can be combined to form circuits. 
Recently, a scalable approach to enzyme-free DNA computing has been pro- 
posed where circuits consist of relatively short DNA strands that communicate 
via strand displacement |31|25| . 

The Poland-Scheraga (PS) model has been very successful in predicting ther- 
mal melting and renaturing of long DNA strands |24|13j . It describes a DNA 



double strand as a ID lattice where each base-pair is either hybridized or open. 
To each state is associated a free energy that has a sequence specific contribu- 
tion from nearest neighbor interactions [29J as well as a polymer contribution 
from the conformational entropy of internal bubbles and frayed ends. Gener- 
alizations of the PS model exists, where the single strands are represented as 
semi- flexible polymers on a 3D lattice |12|14| . This provides a conceptual sim- 
plification since the polymer free energy contributions are given implicitly. The 
Dauxois-Peyrard-Bishop[22j (DPB) model represents a DNA double strand as a 
ID lattice, but each base-pair is described by a continuous base-pair extension. 
The DPB model is defined by a Hamiltonian which includes a hybridization 
potential and a harmonic term penalizing deviations between nearest neighbor 
extensions. 

The chemical structure of short DNA oligomers can be studied with atomistic 
molecular dynamics simulations such as Amber [6|7] and Charmm[4 17J. How- 
ever, if we are interested in mesoscopic properties of long DNA molecules, it 
is more effective to utilize coarse-grained simulation models. Coarse-graining is 
the statistical mechanical process by which microscopic details are systematically 
removed, producing an effective mesoscopic model |16|21j . The major computa- 
tional advantage of coarse-graining is that it allows us to focus our computational 
resources on studying the structures and dynamics at the mesoscopic level. 

Coarse-grained models describe a nucleotide by a small number of effective 
interaction sites. In the "three site per nucleotide" model of de Pablo and co- 
workers, three sites represent the phosphate backbone site, the sugar group, and 
the base, respectively [28 27J. There is also a number of "two site per nucleotide" 
models, e.g. the model of Ouldridge and co-workers |19)20| . where one site rep- 
resents the base and another site the backbone and the sugar ring. Savelyev and 
Papoian [30J have formulated a "one site per nucleotide" model. As the number of 
interaction sites per nucleotide is reduced, the chemical structure is progressively 
lost. In simulations of DNA tagged nanoparticles, even more coarse-grained mod- 
els are used. DNA molecules have been modeled e.g. as semi-flexible polymers 
with attractive sites on each monomer or as a single sticky site that can 
be hybridized with free complementary free sticky sites [18J. While the chemical 
structure of DNA has been completely eliminated, these models still retain the 
DNA sequence specific hybridization effects on nanoparticle self-assembly. 

We are interested in studying the statistical mechanics of hybridizing DNA 
strands and in particular the kinetics of DNA self-assembly and DNA computa- 
tion using a DNA model that is as coarse-grained as possible. We have imple- 
mented a general framework allowing directional bonds to be reversibly formed 
and broken during molecular dynamics simulations [33]. Along with the bonds, 
the angular and dihedral interactions required to model the residual effects of 
chemical structure are also dynamically introduced and removed as dictated by 
the bond dynamics. This framework allows us to simulate reversible hybridiza- 
tion of complementary beads and chains built from such beads. In the present 
paper, we study a minimal dynamic bonding DNA model. For simplicity, we 
assume that the binding energy, as well as the bond, angular, and dihedral po- 



tentials are independent of sequence, and we have chosen a force field that pro- 
duces a flat ladder-like structure in the double stranded state. Our motivation 
for these choices are to minimize the number of parameters required to specify 
the DNA model. 

Dynamic bonding DNA models combine ideas from most of the existing DNA 
models. We regard them as dynamic generalizations of statistical mechanical 
theories and simultaneously as simplifications of coarse-grained DNA models. As 
in the PS model, a complementary base-pair can either be hybridized or open. 
When a base-pair is hybridized, it is characterized by a continuous hybridization 
potential as in the DPB model. Dynamic bonding DNA models can also be 
regarded as off-lattice generalizations of the lattice PS model [12]. Rather than 
trying to model chemical structure with interaction sites as in the "two and 
three site per nucleotide" models |28|27|1 9 20j dynamic bonding DNA models use 
angular and dihedral interactions to model the residual effects of local chemical 
structure. Dynamic bonded DNA double stands can reversibly melt and reanneal, 
which is not possible with the "one site per nucleotide model" of Savelyev and 
Papoian [30J. Finally, as in the sticky DNA models[18j, a single bead in a dynamic 
bonding DNA model can equally well represent a domain. 

Sect. [2] presents the dynamic bonding DNA model, which is used in Sect.[3]to 
study self-assembly of DNA constructs and DNA-computing constructs. Sect. [4] 
ends the article with a conclusions. 

2 Dynamic bonding DNA model 

In the present dynamic bonding DNA model, single stranded DNA (ssDNA) is 
represented by a string of nucleotide beads connected by stiff springs representing 
directional backbone bonds. Instead of using a four letter alphabet representing 
the ACGT nucleotides, in the present paper we increase the alphabet maximally 
to avoid getting trapped in transiently hybridized states. Physically, this corre- 
sponds to assuming that each bead represents a short sequence of nucleotides 
i.e. a domain, and that two non-complementary beads or domains are unable 
to hybridize. A novel feature of our DNA model is that it involves dynamic 
hybridization bonds, which are introduced or removed between complementary 
interaction sites or beads when they enter or exit the hybridization reaction ra- 
dius. Along with the bonds, we dynamically introduce or remove angular and 
dihedral interactions in the chemical neighborhood of a hybridizing bead pair. 
These interactions are introduced based on the local bond and bead type pat- 
tern, and hence allows us to retain some effects of the local chemical structure in 
coarse-grained models. We utilize bonds carrying directionality to represent the 
3'-5' backbone structure of DNA molecules. This allows us to introduce dihe- 
dral interactions that can distinguish between parallel and anti-parallel strand 
alignments. We have implemented this framework in a modified version of the 
Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [23 33J. 

The DNA model relies on two ingredients, a Langevin dynamic for propagat- 
ing a system in time and space, and a dynamic directional bonding scheme [33J 



that propagates the chemical structure of the system. The force on bead i is 
given by a Langevin equation 

Fi = — Vr./7 — — Hi + £i With U = [/"bond + U ang \ e + ^dihedral + ^pair- 

Here, the first term denotes a conservative force derived from the potential U. 
The second term is a velocity dependent friction, and the third a stochastic driv- 
ing force characterized by (£i(t)£j(t')) = feTm / '(r 'At)5ijS(t — t r ) . The potential 
U comprises four terms representing bond, angular, dihedral, and non-bonded 
pair interactions, respectively. The friction and stochastic driving force implicitly 
represents the effect of a solvent with a specified friction and temperature. The 
Langevin dynamics is integrated using a Velocity Verlet algorithm with a time 
step At = O.OOlrx and r = 2tl using a customized version of LAMMPS [ 23|33| . 

Here and in the rest of the paper we use reduced units defined by the Langevin 
dynamics and DNA model. The unit of energy is e = &#T, where we set Boltz- 
mann's constant ^ to unity. The bead-to-bead distance along a single strand 
defines the unit of length a which correponds to the rise distance of DNA. The 
mass is m = 1 for all beads. A Langevin unit of time is defined as tl = ay/m/e. 
The diffusion coefficient of DNA model strand is D(n) = ksTT / (mn) where n 
is the total number of beads. This an be can be equated with the DNA diffusion 
coefficient of a particular experimental conditions to obtain a time mapping. 
Extrapolating the data in Ref. [34J yields tl ~ 1.6 x 10 -12 s for n = 20. 




Figure 1. Illustrative DNA conformation, a) complementary beads, backbone and hy- 
bridization bonds, b) angular interactions indicated by two lines parallel to the involved 
bonds, c) dihedral interactions indicated by three lines parallel to the involved bonds. 
The figure is explained in the text. 

Fig. [IJi shows complementary nucleotide beads with the same hue but dif- 
ferent levels of color saturation. As a simplification, we allow each bead only to 
hybridize with a single complementary bead. The DNA model has two types of 
bond interactions: permanent backbone bonds (shown green/blue) and dynamic 



hybridization bonds (shown red). Backbone bonds and hybridization bonds are 
characterized by the two potentials: 



^o„d )bb (r) = ^% ((r - r*f - (r c b - r b ) 2 ) 

and 



^bond,hyb(V) 



r ^r^((r-r h ) 2 -(r c h -r h )- 2 ) for r < r\ 
for r > r!?. 



In the simulations, we use Z7 m in 5 bb = lOOe, r$ = 1<j, and = 1.2a, rjj = 2a 
and rj? = 2.2a. Note that ?7bond,hyb(0 < for all distances. When two non- 
hybridized beads of complementary type are within a reaction distance r£ a 
hybridization bond is introduced between them. If they move further apart than 
r£ again, the hybridization bond is broken. The pair-interaction between beads 
is given by a soft repulsive potential, while we use the same potential for angular 
and dihedral interactions. They are given by 



( 7rr 

l + cos(^ 
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where we use A = le and r£ = la in the simulations, and 

U{9] 00, tfmin) = (COS[0 - O ] + 1) , 

Along the backbone of single strands we use a permanent angular interaction 
defined by U(0; 0q = 7r, U m i n = 25e). This determines the persistence length of 
single strands. In Fig. [TJd backbone angular interactions are shown as thick red 
lines around the central bead defining the angle. 

In real DNA molecules, the hydrogen bonds between Watson-Crick comple- 
mentary nucleotides act together with stacking interactions and the phospho- 
rdiester backbone bonds to give rise to a helical equilibrium structure of the 
double strand. In our coarse-grained model, we utilize angular and dihedral in- 
teractions to determine the ladder-like equilibrium structure of our DNA model. 
To control the stiffness of the double strands and to ensure anti-parallel 3'-5' 
alignment of the two single strands, we have assigned directionality to the back- 
bone bonds [33]. This is also necessitated by the fact that the 3' and 5' carbons 
of the nucleotide sugar ring have been merged into the single nucleotide bead. 
Fig. [T^l shows the backbone bonds colored green/blue to indicate the 3' and 5' 
ends, respectively. 

When a hybridization bond is introduced, we also dynamically add angu- 
lar interactions between the hybridization bond and the neighboring backbone 
bonds. These angular interactions are characterized by the potential U (0; &o = 
7r/2, C/mi n5 a), which favors a right angle conformation. When a hybridization bond 
is broken, concomitantly all the associated angular interactions are removed. In 
Fig. [TJd the angular interactions are shown as green lines indicating the angle. 



Besides introducing angular interactions, we also dynamically introduce dihe- 
dral interactions. A dihedral interaction involves four beads connected by three 
bonds, which defines a particular bond pattern, where the bonds can either be a 
hybridization bond, a 3' - 5 ; backbone bond, or a 5' - 3' backbone bond. Three 
bond patterns are possible. The bond pattern corresponding to red dihedrals 
in Fig. [I]:, is characterized by U(0; Go = 0, Z7 m i n ,d) which favors a planar (cis) 
conformation. The bond pattern corresponding to blue dihedrals is characterized 
by U(G; 0$ = 7r, /7 m in 5 d, a = 0) which favors parallel backbone (trans) conforma- 
tion. The last dihedral pattern corresponding to green dihedrals is characterized 
by U(G;Oq = 0, E/ m in,d) which favors a parallel (cis) conformation. Note that 
without the directional backbone bonds, we would not be able to distinguish 
between these two latter dihedral patterns. 

During a simulation, at each time we introduce a hybridization bond, we also 
introduce up to four angular interactions and up to eight dihedral interactions, 
less if the hybridization bond is at the end of a strand. Let A be the total decrease 
in binding energy when two beads hybridize inside a chain, and we assign one 
third of this energy to bond, angular, and dihedral interactions, respectively. This 
choice does not affect the static properties of the model, which are determined 
by the total energy associated with a conformation, however it does influence 
the dynamic properties. Hence U m i n ^ y b = A/3, U m i nj3i — Aj 12, and Z7 m in 5 d = 
A/18. We define A — lOe as a reference energy. Since only the ratio A/T enters 
the partition function of the model, this effectively fixes the absolute melting 
temperature of the stands. With the present model, the time spend per particle 
per step is approximately 1 x 10 -5 s on a standard PC. 

3 Results 

Three dimensional DNA structures can be built by utilizing the self-assembly 
properties of complementary strands and by linking several stands into a e.g. 
end-linked constructs. In particular, we have designed four constructs each com- 
prising three end-linked 16 bead long strands. By programming the complemen- 
tarity of the strands, we have designed the constructs to self-assemble into a 
tetrahedron |10|9| . We have also programmed the complementarity of 12 DNA 
constructs each comprising 5 end-linked 8 bead long strands. These constructs 
have been designed to self-assemble into an icosahedron[3j. We estimate that the 
melting temperatures are T m (8) ~ 1.3e, and T m (16) ~ 1.6e from a separate set 
of melting simulations (not shown). 

Fig. [2] shows visualizations of the DNA constructs during the self-assembly 
process. Initially the constructs are randomly placed into the simulation box. 
Progressively, complementary strands hybridize with each other, and the con- 
structs form fragments that ultimately yield the designed target structures. The 
time scale of the self-assembly dynamics is determined by the time it takes the 
constructs to diffuse, collide, and hybridize completely. Since we have the sim- 
ulation trajectory, we can also characterize the detailed time dependence of the 
self-assembly dynamics. Fig. [3] shows the fraction of hybridized bonds as a func- 
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Figure 2. Self-assembly of a tetrahedron from four 3-fimctional DNA constructs (top 
row) and an icosahedron from twelve 5-functional DNA constructs (bottom row). Snap- 
shots correspond to times t = IOOOtl, IO.OOOtl, 20.000tl, 50.000tl steps (top row), 
and for t = IOOOtl, 20.000tl, 30.000tl, 60.000tl steps (bottom row). Simulations 
have been performed at T = l.Oe. Note that periodic boundary conditions apply to the 
simulation box. 
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Figure 3. Fraction of hybridized bonds (top) and number of fragments (bottom) vs. 
time and temperature during the self-assembly of a tetrahedron and an icosahedron. 



tion of time. By analyzing the bond structure, we can furthermore study the 
evolution of the number of structural fragments. For the icosahedron, we see a 
slow increase in the hybridized bond fraction towards unity as the structure is 
progressively assembled, while the number of fragments drops simultaneously 
from initially twelve free constructs to one when all constructs form a single 
icosahedron. However, even through we only have a single fragment, it takes fur- 
ther time for the remaning hybridization bonds to be formed. The equilibrium 
hybridization bond fraction does not appear to be reached at the end of the 
simulation at 1 x 10 8 r. For the tetrahedron, we observe a similar increase in the 
fraction of hybridized bonds, however with six distinct steps corresponding to 
the hybridization of each edge of the tetrahedron. 

The self-assembly dynamics is stochastic and depends on initial conditions 
and random diffusive motion. We have run some of the simulations twice to see 
how they approach equilibrium along different trajectories. The equilibrium hy- 
bridization bond fraction appears to have been reached by the tetrahedron self- 
assembly simulations. For tetrahedra, we observe that self-assembly at higher 
temperatures leads to a marked decrease in the average hybridization bond frac- 
tion similar to melting of DNA double strands. At T = 1.8e the temperature is 
above the melting temperature of the DNA constructs, and they only transiently 
hybridize. Since we have a single fragment at equilibrium, the bond reduction 
is most likely due to DNA bubbles. From the data sets we can estimate that 
the melting temperature of the tetrahedron i.e. 0(T m ) = 0.5 is approximately 
T m (Tetrahedron) « 1.5e. 

In the strand displacement approach to DNA computation, individual gates 
consist of one DNA template that is composed of several logical domains. In 
their initial state, all domains but one are hybridized to one or more comple- 
mentary strands and are therefore inert. The only exposed single strand domain 
of each gate is a short toehold region at one end of the template. This toehold 
region can reversibly bind a complementary signal strand which is designed to 
be longer than the toehold domain and complementary to the next domain (s) 
of the template. The newly binding signal is then able to hybridize to all match- 
ing domains of the template, thereby displacing strands that where previously 
bound [37J. The displaced strands can be nuorescently marked output signals, or 
internal signals that can bind to toehold regions of downstream gates. By choos- 
ing domains of appropriate length, it can be guaranteed that toehold binding is 
reversible, whereas the final strand displacement is effectively irreversible, thus 
computation is energetically downhill and kinetically irreversible, if and only if 
the correct input strands are present and match the logical setup of the gates. 
It has been shown that this approach leads to modular logic gates that enable 
the design of large scale DNA circuits |5|15| . 

Fig. [4] shows simulations of the strand displacement process underlying Seelig 
et al.'s DNA computing approach |31| . The top row shows the successful displace- 
ment of an initially hybridized 12 bead long signal strand from a 20 bead long 
template by a 20 bead long signal strand: once the signal strand diffuses to and 
binds to the toehold region, branch migration occurs quickly (during 300 time 
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Figure 4. Simulations of strand displacement. A 20 bead long oligomer displaces 
a 12 bead long oligomer initially hybridized to a template for times t =500, 1.600, 
1.700, 1.900tl (top row). A 10 + 10 bead long oligomer where the latter half is non- 
complementary fails to displace a short oligomer hybridized to a template for times 
t = 100, 3.000, 7.000, 10.000tl. Simulations are run at temperature T = le. Non- 
complementary beads are show as gray. 
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Figure 5. Time evolution of branch migration on individual template nucleotide beads. 
From top to bottom a) A 12 bead long oligomer (red) being displaced by a 20 bead 
long oligomer (green) on the template, b) a 20 bead long oligomer with different ran- 
dom conditions, c) 12 bead oligomer (red) competing with a 12 + 8 oligomer with 12 
complementary and 8 non-complementary beads (green), d) a 12 bead oligomer (red) 
competing with a 10+10 oligomer, e) a 12 bead oligomer (red) competing with a 8 + 10 
oligomer, f) a 12 bead oligomer (red) competing with 5+10 oligomer. Also shown is 
the branch migration point (black). 



units) and the formerly bound signal strand is displaced irreversibly. The bot- 
tom row, on the other hand, shows how the displacement stalls in the presence of 
mismatches: here, a mismatch in the domain (last 10 beads) permits further hy- 
bridization of the signal strand. The newly binding and the original signal strand 
compete for matching bases in a random walk process until the non-matching 
strand dehybridizes again and leaves the gate available for potential matching 
signals (not present in the simulation). 

Fig. [5] shows statistics of the displacement processes for several runs: the 
graphs depict hybridized bases of the original (red) and the newly binding sig- 
nal strand (green), as well as the branch migration point (black). In the case of 
matching signals (top two simulations), it can be seen that displacement occurs 
quickly and essentially irreversible once the original strand is fully displaced. 
In the third simulation the signal strand and hybridized strand has the same 
length, and the interface is seen to diffuse forwards and backwards. A single 
dehybridization event is also observed for the original strand. In the case of mis- 
matching signals (bottom three simulations), the displacement cannot proceed 
further than nucleotide 10, and the interface randomly moves between posi- 
tions 8 and 10, until - occasionally - the mismatching signal dehybridizes from 
the toehold region (lack of green markers). In this case, the number of beads 
complementary to the toehold region (here 10, 8, and 5 beads) determines the 
equilibrium between hybridized and dehybridized configurations, and thus the 
performance and availability of the gate. Fig. [5] also depicts a source of potential 
failure in logical gates based on strand displacement, as the output signal can 
spontaneously dehybridize even in the absence of a matching input signal (as 
observed in the fourth simulation). 

4 Conclusions 

With these initial simulations, we have demonstrated that our coarse-grained 
DNA model can succesfully simulate DNA assembly as well as DNA strand 
displacement dynamics which form the basis of state-of-the-art DNA computing 
approaches. We have successfully simulated self-assembly of DNA tetrahedra and 
icosahedra from four and twelve branched DNA constructs, respectively. Simula- 
tions show that the constructs self-assemble into the expected target structures. 

We have further simulated successful displacement of an output strand when 
a matching input strand is present. In the presence of mismatches, we could 
demonstrate how the displacement process is prevented. Our simulations also 
capture potential failures of gates based on strand displacement, namely spon- 
taneous release of the output strand in the absence of an input signal. These 
proof of concept simulations demonstrate how our coarse-grained model can be 
used to optimize the length and arrangement of toehold and domain structures 
in DNA computing approaches. 

While such gate optimizations do not necessarily require spatially resolved 
models, our coarse-grained DNA model enables us to study systems that inte- 
grate DNA assembly and computing within a single framework. This enables us 



to use these simulations as a starting point for building and testing statistical 
mechanical theories describing these complex systems. 
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